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Abstract 

Heavy metals released by anthropogenic activities such as mining trigger pro- 
found changes to bacterial communities. In this study we used 16S SSU rRNA 
gene high-throughput sequencing to characterize the impact of a polymetallic 
perturbation and other environmental parameters on taxonomic networks within 
five lacustrine bacterial communities from sites located near Rouyn-Noranda, 
Quebec, Canada. The results showed that community equilibrium was disturbed 
in terms of both diversity and structure. Moreover, heavy metals, especially cad- 
mium combined with water acidity, induced parallel changes among sites via the 
selection of resistant OTUs (Operational Taxonomic Unit) and taxonomic domi- 
nance perturbations favoring the Alphaproteobacteria. Furthermore, under a simi- 
lar selective pressure, covariation trends between phyla revealed conservation and 
parallelism within interphylum interactions. Our study sheds light on the impor- 
tance of analyzing communities not only from a phylogenetic perspective but also 
including a quantitative approach to provide significant insights into the evolu- 
tionary forces that shape the dynamic of the taxonomic interaction networks in 
bacterial communities. 



Introduction 

Anthropogenic activities such as mining and smelting trig- 
ger profound damage to aquatic ecosystems (van Dam 
et al. 2002; Eisler 2004). Rational rehabilitation strategies 
are an urgent requirement to help such ecosystems to 
recover their original functions. However, efforts to restore 
aquatic biomes frequently fail to fully re-establish ecosys- 
tem services (Peralta et al. 2010). Bacterial communities 
are integral service providers (Gutknecht et al. 2006), thus 
conservation and rehabilitation of an impacted area require 
a thorough understanding of how multiple interactive spe- 
cies adapt to face such dramatic environmental perturba- 
tions. 

Metabolic plasticity is likely to underpin the fundamental 
contribution bacterioplankton make to aquatic ecosystems 
(Armitage et al. 2003; Comte and del Giorgio 2011). As 
such, environmental fluctuations elicit a rapid adaptive 
response at the community level. However, the processes 
underlying this metabolic plasticity remain unclear. There 



is some evidence that metabolic pathways of lacustrine bac- 
terioplankton may simply be compartmentalized between 
the main phyla (Actinobacteria, Betaproteobacteria, Alpha- 
proteobacteria, and Bacteroidetes) (Debroas et al. 2009). 
Such functional compartmentalization suggests that inter- 
acting networks of taxa may constrain the functional reper- 
toire of the whole microbial community. 

Where taxonomic functional compartmentalization 
exists, so does a suite of mechanisms to mitigate the impact 
of the loss of species richness from the microbial commu- 
nity. Horizontal gene transfer (HGT) is universally 
reported as the major mechanism to facilitate rapid adapta- 
tion in microbial communities (Shintani et al. 2008; So- 
becky and Hazen 2009; Parnell et al. 2010; Zhaxybayeva 
and Doolittle 2011). The resultant fluidity of functional 
characteristics between taxa can potentially bypass keystone 
species within functional networks. However, HGT fre- 
quency and efficiency is known decrease with genetic diver- 
gence (Andam and Gogarten 2011) and this phenomenon 
would in turn limit community level adaptation where 
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functional characteristics are discretely packaged among 
only distantly related taxa. 

Moderate and transient environmental changes do not 
dramatically reduce the catalog of ecosystem services pro- 
vided by bacterial communities even though metabolic 
function maybe compartmentalized among phyla (de Bello 
et al. 2010). In fact, at intermediate levels of disturbance, 
diversity is maximized because both competitive K-selected 
and opportunistic R-selected species can coexist (Dial and 
Roughgarden 1988; Kassen et al. 2004). Alternative mem- 
bers of the community with similar functional characteris- 
tics, replace the lost species (Folke et al. 1996; Lee et al. 
2010; Barberan et al. 2012). As such, a reticulate network 
of community-level functional interactions buffers moder- 
ate environmental perturbations (Laplante and Derome 
2011). In contrast, strong and permanent environmental 
disturbances are expected to dramatically alter bacterial 
community diversity (Rohr et al. 2006; Deng et al. 2007; 
Parnell et al. 2009). Such perturbations, for example the 
presence of an industrial pollutant, are expected to erase 
many species either via the direct action of the pollutant, or 
- crucially - via secondary effects e.g. the loss of keystone 
species or phyla from the functional network. 

Mining and smelting release toxic heavy metals through 
acid mine drainage (AMD). AMD severely impacted sites 
are extreme environments as very acidic, very low in 
organic matter and rich in heavy metals or metalloids. 
When in excess in living cells, trace metals can bind to 
physiologically sensitive intracellular molecules or organ- 
elles, leading to deleterious effect (Wallace et al. 2003). 
Polymetallic contamination restricts community member- 
ship to organisms well-adapted to strong acidity and ele- 
vated concentration of metal sulfides (Collon 2003; Shade 
and Handelsman 2011). In this respect, heavy metals exert 
a very strong selective pressure on bacterial communities, 
causing structural changes in the community itself and are 
usually accompanied with a reduction in the microbial bio- 
mass (Kandeler et al. 2000), and with a loss of biodiversity 
(Deng et al. 2007; Laplante and Derome 2011). Moreover, 
communities exposed to environmental stress like radionu- 
clide or polymetallic contamination typically share more 
structural similarities when compared with communities 
from unpolluted lakes (Fields et al. 2005; Laplante and 
Derome 2011). Studying these patterns of convergent or 
parallel evolution provides an opportunity to identify key 
associated parameters (Schluter 2000; Derome and Bernat- 
chez 2006; Derome et al. 2006). In a previous survey we 
detected clear evidence of parallel adaptation in a lake sys- 
tem of both connected and independent bacterioplankton 
communities using 16S-SSU-rDNA DGGE profiling (Lap- 
lante and Derome 2011). Initially, we showed that polyme- 
tallic contamination drove parallel adaptation at the 
taxonomic level in two independent bacterial communities. 
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Furthermore, we also showed that long-term exposure to 
contaminant gradients drove significant taxonomic diver- 
sity changes within three interconnected bacterial commu- 
nities (Laplante and Derome 2011). However, 16S-SSU- 
rDNA DGGE fingerprint techniques are insufficiently sen- 
sitive to allow quantification of species abundance and pro- 
vide no access to the rare biosphere. In contrast, next 
generation amplicon deep sequencing of the 16S locus pro- 
vides much greater sensitivity. Taxonomic diversity and 
taxon abundance can be simultaneously evaluated with 
high reproducibility (Pilloni et al. 2012). Such accurate 
microbial community profiles enable downstream network 
and systems analysis of interactions between taxa. 

In the following experiment, barcoded 454-pyrosequenc- 
ing was used to study the same lake system as in Laplante 
and Derome (2011), to provide deeper insight into (i) the 
membership and structural changes occurring in the bacte- 
rial communities under long-term exposure to polymetallic 
contaminants, (ii) the parallel changes occurring between 
two independent bacterial communities impacted by simi- 
lar selective pressure and (iii) the related environmental 
process driving community changes in the studied system. 

Material and methods 

Study area and sampling sites 

Detailed study area and sampling site information was 
detailed in Laplante and Derome (2011). Briefly, five sam- 
pling sites were chosen in the surroundings of Rouyn- 
Noranda (Quebec, Canada) (see Fig. 1), where a strong 
mining activity generated polymetallic contaminant gradi- 
ents (Couillard et al. 1993; Giguere et al. 2003, 2006). 
Opasatica Lake (Opa) was chosen as an unpolluted lake 
reference, while Turcotte Lake (Tur) was chosen as a con- 
taminated reference. Hydrologically interconnected sites 
studied were Dasserat Lake (Das), Arnoux Bay (Bar) and 
Arnoux Lake (Lar). Through their interconnection, 
the water flow spreads the AMD from Arnoux Lake to 
Dasserat Lake such that a polymetallic gradient is naturally 
generated. 

Sampling and filtration procedures 

These procedures were described in Laplante and Derome 
(2011). Briefly, sampling was done on June 2010 by collect- 
ing 6 L of water (3 L in duplicates) at 60 cm of depth in 
the water column. Water samples were filtered first with a 
3.0 /mi mesh size, followed by a 0.22 ^m nitrocellulose 
membrane (Advantec) using a peristaltic filtration equip- 
ment (Masterflex L/S Pump System with Easy-Load II 
Pump Head; Cole-Parmer, Vernon Hills, IL, USA). Dupli- 
cates of the 0.22 /im membranes from the sampled lakes 
were placed into cryotubes containing 1 ml of sterile lysis 
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55°C over night with agitation. All the aqueous phase was 
transferred into a clean microcentrifuge tube containing 
600 ji~L of 6m NaCl (final concentration 2.25 m), mixed and 
centrifuged 20 min at 17 000 g. The supernatant was trans- 
ferred again into a clean microcentrifuge tube containing 1 
volume of ice-cold isopropanol, mixed gently and stored 
30 min at — 20°C. The mixture was centrifuged 20 min at 
16 162 g and the supernatant was thrown away. The pellet 
was washed with ice-cold 70% ethanol, air-dried and finally 
resuspended in 25 /jL of sterile MilliQ H 2 0. Subsequently, 
DNA purity (260/280 ratio) and quantity were measured 
using a Nanodrop instrument (ND-1000, Nanodrop, 
Thermo Fisher Scientific, Waltham, MA, USA). 






c 


Home smelter 




0 


Aldermac mine site 




■ 


Rouyn Noranda 




■ 


Sampling site 


0 




13km 



Figure 1 Geographical localization of Rouyn-Noranda (Abitibi-Temisc- 
amingue, Canada) and the sampling sites visited in June 201 0. Contam- 
inated reference site: Turcotte Lake (Tur); clean reference site: 
Opasatica Lake (Opa); tested lake system: Arnoux Lake (Lar), Arnoux 
Bay (Bar), and Dasserat Lake (Das). 



buffer (40 mM EDTA, 50 mM Tris-HCl, 0.75 m sucrose) 
and then flash frozen into liquid nitrogen until DNA 
extraction. 



DNA extraction 

Genomic DNA was extracted using a modified protocol of 
salt-extraction from Aljanabi and Martinez (1997). During 
the first lysis step, 113 of lysozyme (final concentration 
1 mg/mL) was added to the sample and incubated for 
45 min at 37°C. After this step, 113 of proteinase K (final 
concentration 0.2 mg/mL) and 450 jiL of SDS 10% (final 
concentration 1%) were added to the lysate and incubated at 



Barcoded 16S pyrosequencing 

Each sample, consisting of 12 ng of extracted genomic 
DNA, was PCR amplified using Takara Extaq Premix 
(Fisher). All PCR reactions were performed in a final vol- 
ume of 50 /iL containing 25 fiL of Premix Taq, 1 of 
each primer and sterile MilliQ H 2 0 to up to 50 /(L. To 
achieve the PCR amplifications, a general reverse primer 
(R519) combined with B primer (Roche) was used in com- 
bination with a unique tagged forward primer (F63-tar- 
geted) combined with A primer (Roche) (see Table 1). 
PCR conditions were as follow: after a denaturing step of 
30 s at 98°C, samples were processed through 30 cycles 
consisting of 10 s at 98°C, 30 s at 55°C, and 30 s at 72°C. 
The final extension step was done at 72°C for 4 min 30 s. 
Following amplification, samples were purified using 
AMPure Beads (Beckman Coulter Genomics, Denver, MA, 
USA). Samples were adjusted to 100 fiL with EB (Qiagen, 
Hilden, Germany). Then, 63 /iL of beads were added, the 
samples were mixed and incubated for 5 min at room tem- 
perature. Using a Magnetic Particle Concentrator (MPC), 
the beads were pelleted against the wall of the tube and the 
supernatants were removed. The beads were washed twice 
with 500 ,uL of 70% ethanol and incubated for 30 s each 
time. The supernatants were removed and the beads were 
allowed to air dry for 5 min. The tubes were removed from 
the MPC and 24 iiL of EB were added. The samples were 
vortexed to resuspend the beads. Finally, using the MPC, 



Table 1 . Details of the primers used for the 1 6S barcoded pyrosequencing. 



Primer 


Adaptor sequence 


Barcode sequence 


Primer sequence 


R519 

F63_Lar 

F63_Bar 

F63_Das 

F63_Tur 

F63_Opa 


CCTATCCCCTGTGTGCCTTGGCAGTCT 

CCATCTCATCCCTGCGTGTCTCCGACTCAG 

CCATCTCATCCCTGCGTGTCTCCGACTCAG 

CCATCTCATCCCTGCGTGTCTCCGACTCAG 

CCATCTCATCCCTGCGTGTCTCCGACTCAG 

CCATCTCATCCCTGCGTGTCTCCGACTCAG 


ACGAGTGCGT 
ACGCTCGACA 
AGACGCACTC 
AGCACTGTAG 
ATCAGACACG 


CAGGWATTACCGCGGCKGCTG 

CAGGCCTAACACATGCAAGTC 

CAGGCCTAACACATGCAAGTC 

CAGGCCTAACACATGCAAGTC 

CAGGCCTAACACATGCAAGTC 

CAGGCCTAACACATGCAAGTC 


Lar, Arnoux Lake; Bar, Arnoux Bay; Das, Dasserat Lake; Tur, Turcotte Lake; Opa, Opasatica Lake. 
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the beads were pelleted against the wall once more and the 
supernatants were transferred to a new clean tube. The 
samples were quantified with Nanodrop and mixed equally 
then sent to the Plateforme d' Analyses Biomoleculaires 
(Institut de Biologie Integrative et des Systemes, Universite 
Laval) to perform the pyrosequencing using the 454 GS- 
FLX DNA Sequencer with the Titanium Chemistry (Roche, 
Basel, Switzerland) using the procedure described by the 
manufacturer. 

Pyrosequencing data analysis 

The open-source community-supported software program, 
Mothur (http://www.mothur.org) (Schloss et al. 2009), 
was used to process and analyze the sequence data follow- 
ing the Costello Stool Analysis tutorial (http://www. 
mothur.org/wiki/Costello_stool_analysis) (Costello et al. 
2009). Following this tutorial, the filtration, the alignment 
and the trimming of the sequences was done to extract the 
high-quality reads at a uniform length of 300 nt to facilitate 
comparative analyses. Accordingly, the data were filtered 
using the following criteria: minimum average quality 
score: 35 for a windows of 50, number of differences to the 
primer sequence = 0, maximum number of differences to 
the barcode sequence = 0, number of ambiguous base 
calls = 0, maximium homopolymer length = 8. Then 
trimmed reads were grouped into clusters with a 97% iden- 
tity threshold. The Costello Stool Analysis tutorial also pro- 
vided the commands required to (i) calculate diversity and 
richness indexes (Nonparametric Shannon Index, Simpson 
Index, sequences coverage, abundance of unique OTUs, 
Chaol estimate of richness), (ii) achieve the taxonomical 
survey of each community, (iii) build Venn diagrams at a 
cutoff of 3%, (iv) construct dendrograms showing either 
the taxonomical diversity similarities between communities 
(using the laccard Coefficient) or their taxonomical struc- 
ture (i.e. OTU relative abundances) (using the Theta Coef- 
ficient) (Yue and Clayton 2005), and (v) determine 
whether the clustering within the trees is statistically signifi- 
cant using the Unifrac unweighted and weighted methods 
(Lozupone and Knight 2005). 

Statistical analysis 

All data were first tested for normality and equality of vari- 
ances (Zar 1999). Linear regressions were performed using 
the R command 'lm(tableau$X~tableau$Y)' to test the exis- 
tence of a linear relation between the response variable Y, 
namely the diversity indices and the taxonomical surveys of 
each community, and the explicative variable X, namely the 
abiotic data collected on the field (trace metals [Al, Cd, Cu, 
Fe, Mn, Pb, Zn], major cations [Ca, Mg, Na, K, S], DOC, 
pH and temperature). Then, analyses of variances (anovas) 
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were undertaken using the R commands 'model = lm(tab- 
leau$X~tableau$Y)' and 'anova(model)' to assess the fac- 
tors explaining the variance of the models. Effects for the 
anova tests were deemed significant at P < 0.05 and mar- 
ginally significant at P = 0.05—0.10. In addition, for each 
linear regression, the coefficient of determination (R 2 ) was 
calculated to assess the validity of the model. The statistical 
analyses were performed using R statistical software (http:/ 
www.r-project.org/). 

Bacterial network analysis 

Based on the relative abundance of the 50 dominant spe- 
cies, the concentrations of cadmium, pH, and dissolved 
organic carbon (DOC) values measured in the five lakes, 
bacterial networks were built using Rcmdr, a platform- 
independent menu interface running with R software. The 
correlation coefficients were calculated by applying the 
Spearman algorithm. This previous filtering step removed 
poorly represented OTUs and reduced network complex- 
ity. We considered a valid co-occurrence event to be a 
robust correlation if the Spearman's correlation coefficient 
was both >0.6 and statistically significant (P-value <0.01) 
(Barberan et al. 2012). The networks were displayed using 
the Fruchterman-Reingold algorithm with 500 iterations. 
The nodes in the reconstructed networks represent the 
OTUs with 97% identity and the selected physicochemical 
parameters. Then, the edges (that is, connections) corre- 
spond to a strongly significant interaction between nodes 
and red edges indicate a negative interaction between two 
nodes. 

Results 

Environmental data 

Trace metals (Al, Cd, Cu, Fe, Mn, Pb, Zn), major cations 
(Ca, Mg, Na, K, S), DOC, pH, and temperature for the 
studied lakes are summarized in Table 2. For more details 
about the abiotic survey procedure, see Laplante and Der- 
ome (2011). 

454 data analyses 

After trimming the primer sequences, barcodes and adapter 
tags, the average sequence length was 425 nt and the total 
dataset comprised 47 124 high quality sequences with a 
minimum length of 300 nt. Then, after removing the chi- 
meras and retrieving the unique sequences, the analyses 
using Mothur generated multiple outputs, which are sum- 
marized in Table 3. The Nonparametric Shannon Index, 
which takes into account the fact that possible rare species 
were not collected during the sampling, was the highest in 
Opa (H' = 4.564), followed by Das (H' = 3.834), Tur 

© 2013 The Authors. Published by Blackwell Publishing Ltd 6 (2013) 643-659 
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Table 2. Abiotic parameters measured at each sampling site. 

Detection limit Lar Bar Das Tur Opa Mean ± SD 

Trace metals (mg/L) 

Fe 0.002 4.939 0.705 0.11 0.097 0.07 1.184 ± 2.116 

Al 0.001 0.930 0.405 0.052 0.113 0.074 0.315 ±0.372 

Zn 0.0007 0.5574 0.2703 0.0347 0.1063 O.0007 0.1939 ± 0.2282 

Mn 0.0001 0.4421 0.2598 0.0036 0.0368 0.0019 0.1488 ± 0.1960 

CU 0.0005 0.0507 0.0243 0.0075 0.0141 0.0029 0.01 99 ± 0.01 90 

Pb 0.003 0.003 <0.003 <0.003 <0.003 <0.003 0.003 ± 7.91 6 E-05 

Cd 0.0002 0.0010 0.0007 0.0002 0.0009 O.0002 0.0006 ± 0.0004 

Major cations (mg/L) 

S 0.02 22.03 13.16 5.10 2.33 1.93 8.91 ± 8.61 

Ca 0.02 9.92 8.26 7.34 2.00 7.97 7.10 ± 3.00 

Mg 0.002 3.490 2.675 1.987 0.434 2.452 2.207 ±1.131 

Na 0.01 1.26 1.24 1.21 0.68 3.24 1.53 ±0.99 

K 0.002 0.597 0.461 0.497 0.139 0.917 0.522 ± 0.280 

Others 

Doc 0.5 4.4 1.8 7.2 3.8 7.7 5.0 ± 2.5 

pH 0.05 3.77 4.69 7.11 4.91 7.64 5.62 ±1.67 

Temp. 0.5 19.0 17.5 17.0 17.0 16.5 17.4 ±1.0 



Lar, Arnoux Lake; Bar, Arnoux Bay; Das, Dasserat Lake; Tur, Turcotte Lake; Opa, Opasatica Lake; SD, Standard Deviation; Al, Aluminum; Ca, Calcium; 
Cd, Cadmium; Cu, Copper; Fe, Iron; K, Potassium; Mg, Magnesium; Mn, Manganese; Na, Sodium; Pb, Lead; S, Sulfur; Zn, Zinc; DOC, Dissolved 
Organic Carbon; Temp., Temperature. 



Table 3. Data collected for each community using Mothur. 



Site A/ seqs Coverage (%) Shannon InvSimpson S obs S chl 



Lar 12 026 96.7 2.161 0.713 532 1794.75 

Bar 9819 98.0 0.728 0.141 280 836.32 

Das 8503 92.8 3.834 0.896 847 2674.01 

Tur 5898 94.7 2.924 0.805 444 1204.24 

Opa 1721 82.9 4.564 0.936 405 1308.44 



Lar, Arnoux Lake; Bar, Arnoux Bay; Das, Dasserat Lake; Tur, Turcotte Lake; Opa, Opasatica Lake; W seqs Number of Sequences; Shannon, Nonparamet- 
ric Shannon Index; InvSimpson, Inversed Simpson Index; S 0 t, s , Number of OTUs observed and S C hao. Richness index. 



(H' = 2.924), and Lar (JT = 2.161), while the lowest Shan- 
non value was calculated for Bar (H' = 0.728). The esti- 
mated species richness, as seen with the Chao values, was 
the highest in Das (Schao = 2674.010), followed by Lar 
(Schao = 1794.754), Opa (S ch ao = 1308.438), and Tur 
(Schao = 1204.238), while the lowest richness was observed 
inBar(S C hao = 836.324). 

Relationship between lacustrine microbial communities 

Venn diagrams of the percentages of species shared 
between communities (see Fig. 2) reveal that 4.56% of the 
OTUs were common between Lar and Bar, 2.48% were 
common to Bar and Das, 0.65% were shared between Das 
and Lar, while 0.78% of the OTUs were common to these 
three interconnected sites (see Fig. 2A). By analyzing the 
interconnected sites system with the reference communities 
(see Fig. 2B and C), we found that the unpolluted lakes 



Das and Opa exclusively shared the highest number of 
OTUs with 5.51%, followed by the moderately disturbed 
site Bar and the highly contaminated site Lar with 4.43% 
and 3.45% of exclusively shared OTUs when compared 
with either Opa or Tur, respectively. In contrast, Lar and 
Opa showed the lowest percentage of exclusively shared 
OTUs with 0%. 

The membership and structural similarities between the 
studied communities are illustrated by the distinct dendro- 
grams provided using Mothur (see Fig. 3). Firstly, the 
membership dendrogram (see Fig. 3A), obtained from the 
unweighted (equally weighted) similarity matrix, showed 
two differentiated clusters with a Unifrac value at the node 
of 0.9455. The first cluster clearly regrouped the disturbed 
communities (DC) of Tur, Lar, and Bar (value at the node 
of 0.9313 and 0.0095 substitutions per alignment posi- 
tions), while the second cluster strongly identified the 
undisturbed communities (UC) of Das and Opa (value at 
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Figure 2 Venn diagrams of the percentage of species shared between communities at distance 0.03. The Venn diagrams were obtain using Mothur 
and allowed to primarily assess the percentage of species shared within (A) the tested lake system consisting in Arnoux Lake, high metallic contami- 
nated communities group (Lar, HMC), Arnoux Bay intermediate metallic contaminated communities group (Bar, IMC) and Dasserat Lake, low metallic 
contaminated communities group (Das, LMC), and secondly to assess the percentage of species shared between the tested lake system and (B) the 
clean reference site Opasatica Lake (Opa, LMC) or (C) the contaminated reference site Turcotte Lake (Tur, IMC). 
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Figure 3 UPGMA dendrograms representing the membership and the 
composition shared between communities. The UPGMA dendrograms 
allowed to assess a) the membership (taxonomic diversity) shared 
between communities based on Jaccard coefficient (jclass) and b) the 
composition (taxonomic structure) shared between communities based 
on Yue & Clayton theta coefficient (thetayc). Branch lengths are dis- 
played in black and values at the nodes are displayed in red. In (A), the 
communities clustered into two groups, either the disturbed communi- 
ties (DC) group or the undisturbed communities (UC) group. In (B), fol- 
lowing the polymetallic contamination magnitude over the taxonomic 
structure, the communities clustered into three groups, either the low 
metallic contaminated (LMC) communities group, the intermediate 
metallic contaminated (IMC) communities group or the high metallic 
contaminated (HMC) communities group. Contaminated reference site: 
Turcotte Lake (Tur); clean reference site: Opasatica Lake (Opa); tested 
lake system: Arnoux Lake (Lar), Arnoux Bay (Bar) and Dasserat Lake 
(Das). The scale bar represents 0. 1 substitutions per alignment position. 

the node of 0.8936 and 0.0276 substitutions per alignment 
positions). Secondly, by evaluating the abundance of the 
different taxa, the weighted structural dendrogram (see 
Fig. 3B) indicated three differentiated clusters consistent 
with the magnitude of the contamination. The first cluster 
was the highly metallic contaminated (HMC) community 
of Lar (see Table 2 for detailed abiotic survey), which 



showed a significant unique taxonomical structure (value 
at the node of 0.4723 and 0.4385 substitutions per align- 
ment positions). The other communities were grouped 
together under a similar cluster (value at the node of 
0.4778 and 0.0458 substitutions per alignment positions) 
which comprised Tur and Bar under a subcluster, thus 
being identified as intermediate metallic contaminated 
(IMC) communities (value at the node of 0.5643 and 
0.1011 substitutions per alignment positions from each 
other), and the low metallic contaminated (LMC) commu- 
nities of Opa and Das under a highly significant subcluster 
(value at the node of 0.2916 and 0.2262 substitutions per 
alignment positions from each other). Evident from the 
nodes of both dendrograms, and despite the geographical 
connections between Das, Lar, and Bar, the Das commu- 
nity is more similar to the independent and undisturbed 
reference community at Opa. 

Principal coordinate analyses (PCoA) of the unweighted 
(diversity) and weighted (structure) UniFrac distance 
matrix are showed in Fig. 4. Considering the diversity level, 
PCoA of the unweighted Unifrac distance matrix provided 
compelling evidences about the clustering of the undis- 
turbed communities (UC), namely Das and Opa, and the 
clustering of the disturbed communities (DC) Lar-Bar-Tur. 
Through PCoA, the weighted Unifrac distance matrix 
clearly demonstrated that the communities' structure clus- 
tered in three groups consistent with the magnitude of the 
metallic contamination. The first node clustered the LMC 
communities Opa and Das, while the second cluster closely 
grouped the IMC communities Tur and Bar. Considering 
its unique taxonomical structure, Lar was the sole represen- 
tative of the third cluster identified as a HMC community. 

Taxonomical survey at the diversity level 

Based on each phylum's relative OTUs abundance, the 
taxonomical survey at the diversity level (see Fig. 5) 
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Figure 4 Communities clustered at the membership level and the com- 
position level using PCoA of the unweighted UniFrac distance matrix. 
Communities clustered at (A) the membership level (taxonomic diver- 
sity) and (B) the composition level (taxonomic structure) using PCoA of 
the unweighted UniFrac distance matrix. In (A), two groups are showed, 
the disturbed communities (DC) and the undisturbed communities 
(UC). In (B), three groups are showed following the polymetallic con- 
tamination magnitude, namely the low metallic contaminated (LMC) 
communities, the intermediate metallic contaminated (IMC) communi- 
ties and the high metallic contaminated (HMC) communities. Contami- 
nated reference site: Turcotte Lake (Tur); clean reference site: Opasatica 
Lake (Opa); tested lake system: Arnoux Lake (Lar), Arnoux Bay (Bar), 
and Dasserat Lake (Das). 

revealed that the Proteobacteria were prevalent among the 
disturbed communities. Deeper analysis at the class level 
among the Proteobacteria exerts the high relative propor- 
tion of OTUs belonging to the Alphaproteobacteria 
(Lar = 64.30% > Tur = 56.93% > Bar = 33.58% > Opa = 
30.66% > Das = 15.43%), followed by the Betaproteobacte- 
ria (Bar = 34.34% > Das = 25.76% > Opa =16.25% >Lar 
= 16.17% > Tur = 13.20%). The second most dominant 
phylum based on the relative OTUs abundance was repre- 
sented by the Bacteroidetes, which are found primarily in 
the undisturbed sites and mainly encompassed the classes 
Plavobacteria (Das = 18.35% > Opa = 8.47% > Tur = 
0.65% > Bar = 0.38% > Lar = undetected) and Sphingo- 
bacteria (Das = 9.48% > Opa = 7.78% > Tur = 5.41% > 
Bar = 3.02% > Lar = 0.59%). The third phylum that 
showed high degree of OTU diversity was Actinobacteria 
which, as for the Bacteroidetes, exhibited higher OTU diver- 
sity within the noncontaminated sites. Among the Actino- 




40% 50% 60% 



Figure 5 Variation of the principal phyla relative OTUs abundance 
between communities. Clean reference site: Opasatica Lake (Opa); con- 
taminated reference site: Turcotte Lake (Tur); tested lake system: Ar- 
noux Lake (Lar), Arnoux Bay (Bar), and Dasserat Lake (Das). Phyla 
regrouped under the term Others are the following: Opa = Chlorobi, 
Gemmatimonadetes, Lentisphaerae, Spirochaetes, and Synergistetes; 
Tur = Caldiserica and Chlorobi; Das = Deferribacteres, Gemmatimona- 
detes, and Lentisphaerae; Bar = Aquificae arid Planctomycetes; 
Lar = Chlorobi, Nitrospira, and Spirochaetes. 



bacteria, most representatives were regrouped in the class 
of the Actinomycetales (Das = 16.04% > Opa = 10.30% > 
Bar = 3.40% > Lar = 1.18% > Tur = 0.43%) followed by 
the class Acidimicrobiales (Das = 1.58% > Lar = 0.79% > 
Bar = 0.75% > Tur = 0.43% > Opa = 0.23%). Several 
other phyla were identified within the studied system such 
as the Acidobacteria, Firmicutes, TM7, Chloroflexi, OP10 
and Verrucomicrobia, but they showed very weak relative 
OTU abundances (less than 2.29% in at least one site). 

The relative OTUs abundance survey permitted establish- 
ing the ten most dominant genera at the diversity level 
within each community (see Table 4). The genus Polynucle- 
obacter (class Betaproteobacteria) showed greatest relative 
abundance (Bar =28.30% > Lar = 11.83% > Das = 8.87% 

> Tur = 6.06% > Opa = 5.26%). The second best repre- 
sented genus was Kozakia from the class Alphaproteobacteria 
(Lar = 19.33% > Bar = 0.38% > Das = 0.24% > Tur- 
Opa = undetected), followed by the Acidisphaera 
(Lar = 13.21% > Bar = 1.89% > Tur = 0.87% > Das = 
0.12% > Opa = undetected), the genus Flavobacterium 
from the Bacteroidetes phylum (Opa = 8.47% 

> Das = 6.93% > Bar = 0.38% > Tur-Lar = undetected), 
and the genus Citreimonas from the class Alphaproteo- 
bacteria (Tur = 14.07% > Opa = 0.23% > Das-Bar-Lar = 
undetected). 



Taxonomical survey at the structure level 

Based on each phylum's relative sequences abundance, the 
taxonomical survey at the structure level (see Fig. 6) 
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Table 4. Survey of the ten dominant genera found within each community considering the taxonomical diversity (based on the relative OTUs abun- 
dance). 



Genera 


Relative OTUs 
abundance (%) 


Site 


Genera 


Relative 
abundai 


Kozakia 


19.33 


Tur 


Citreimonas 


14.07 


Acidisphaera 


13.21 




Fodinicurvata 


7.36 


Polyn ucleobacter 


11.83 




Polynucleobacter 


6.06 


Acidocella 


10.06 




Tanticharoenia 


2.81 


Rhodovarius 


3.35 




Sediminibacterium 


1.95 


Tanticharoenia 


2.56 




OP10 Unclassified 


1.73 


Rhodopila 


2.17 




Ponticoccus 


1.52 


Legionella 


1.58 




Wenxinia 


1.52 


Acidobacteria_GP1 


0.79 




Agromonas 


1.30 


Acidiphilium 


0.79 




Rhodopila 


1.08 


Polyn ucleobacter 


28.30 


Opa 


llumatobacter 


5.49 


Novosphingobium 


4.15 




Polynucleobacter 


5.26 


Cohaesibacter 


3.40 




Oryzihumus 


4.12 


TM7 Unclassified 


2.64 




Flavobacterium 


3.20 


Acidocella 


1.89 




Legionella 


3.20 


Acidisphaera 


1.89 




Rhodoferax 


2.52 


Beijerinckia 


1.51 




Novosphingobium 


2.29 


Arsenicicoccus 


1.13 




Pelagibacter 


1.83 


Rhodoferax 


1.13 




Zunongwangia 


1.37 


Nubsella 


1.13 




Laribacter 


1.37 


Polyn ucleobacter 


8.87 








Flavobacterium 


6.93 








Zunongwangia 


4.62 








Arcicella 


3.40 








Lapillicoccus 


3.04 








Persicobacter 


2.79 








Fodinibacter 


2.43 








Rhodoferax 


1.70 








Fulvivirga 


1.70 








Intrasporangium 


1.58 









Lar, Arnoux Lake; Bar, Arnoux Bay; Das, Dasserat Lake; Tur, Turcotte Lake; Opa, Opasatica Lake. 



revealed that the Proteobacteria are dominant among dis- 
turbed communities (Lar = 97.94% > Bar = 98.21% > 
Tur =92.53% > Opa = 70.86% > Das = 61.16%). Deeper 
analysis within the Proteobacteria phylum revealed the high 
relative proportion of sequences belonging to the Betapro- 
teobacteria (Bar = 94.82% > Tur = 36.13% > Das = 

35.06% > Opa = 27.45% > Lar = 19.79%), followed by 
the Alphaproteobacteria (Lar = 77.51% > Tur = 54.68% > 
Opa = 38.99% > Das = 24.02% > Bar = 2.94%). The sec- 
ond most dominant phylum based on the relative 
sequences abundance was represented by the Actinobacte- 
ria, which contained higher OTU diversity within the non- 
contaminated sites (Das = 20.66% > Opa = 16.90% > 
Lar = 1.65% > Bar = 1.37% > Tur = 0.20%). Among the 
Actinobacteria, most representatives were regrouped in the 
class of the Actinomycetales (Das = 19.60% > Opa = 
2.06% >Lar = 1.61% > Bar = 1.33% > Tur = 0.17%). 

The third phylum exposing the broader structural com- 
plexity was the Bacteroidetes, which were found mostly in 



the undisturbed sites (Das = 17.78% > Opa = 9.44% > 
Tur = 3.74% > Bar = 0.16% > Lar = 0.15%) and mainly 
regrouped two classes: Flavobacteria (Das = 12.22% 
> Opa = 4.08% > Tur = 0.05% > Bar = 0.02% > Lar = - 
undetected) and Sphingobacteria (Das = 4.85% > Opa = 
2.97% > Tur = 3.64% > Lar = 0.14% > Bar = 0.13%). 

The survey of the relative sequences abundance allowed 
establishing the ten most structurally dominant genera 
within each community (see Table 5). The genus Polynu- 
cleobacter (class Betaproteobacteria) prevailed within all the 
communities (Bar = 94.55% > Tur = 34.34% > Das = 
28.74% > Lar = 19.55% > Opa = 16.61%). The second 
genus with the highest representatives was Kozakia from 
the class Alphaproteobacteria (Lar = 47.66% > Das = 
3.53% > Bar = 1.02% > Tur-Opa = undetected), followed 
by the genus Pelagibacter (Opa = 22.09% > Das = 
18.57% > Bar = 5.09% > Das-Tur = undetected) and the 
genus Fodinicurvata (Tur = 33.37% > Bar = 1.02% > 
Das-Bar-Lar = undetected). 
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Figure 6 Variation of the principal phyla relative sequences abundance 
between communities. Clean reference site: Opasatica Lake (Opa); con- 
taminated reference site: Turcotte Lake (Tur); tested lake system: Ar- 
noux Lake (Lar), Arnoux Bay (Bar), and Dasserat Lake (Das). Phyla 
regrouped under the term Others are the following: Opa = Chlorobi, 
Gemmatimonadetes, Lentisphaerae, Spirochaetes, and Synergistetes; 
Tur = Caldiserica and Chlorobi; Das = Deferribacteres, Gemmatimona- 
detes, and Lentisphaerae; Bar = Aquificae and Planctomycetes; 
Lar = Chlorobi, Nitrospira, and Spirochaetes. 



Physicochemical correlation analyses at the diversity level 

All the correlation results gained from anova analyses 
between the diversity data collected and the physicochemi- 
cal parameters measured at each site are summarized in 
Tables 3 and 6. Among the diversity indices recovered 
using Mothur, only Shannon indices revealed a significant 
correlation with the dissolved organic carbon (DOC) mea- 
sures (R 2 = 0.90; P = 0.01379). 

At the phylum level, statistical analyses also revealed that 
the relative abundance of OTUs belonging to the Actinobac- 
teria was negatively influenced by the cadmium gradient 
(negative correlation, R 2 = 0.95; P = 0.005365) while 
favored by near neutral pH measures (R 2 = 0.87, 
P = 0.02067). Interestingly, the most diverse order within 
the class Actinobacteria, namely the Actinomycetales, mainly 
explained such a negative correlation with cadmium (nega- 
tive correlation, R 2 = 0.88; P = 0.01725). On the reciprocal, 
the diversity of the phylum Proteobacteria was significantly 
correlated with cadmium concentrations (R = 0.90; 
P = 0.01417) and negatively influenced by alkaline pH envi- 
ronments (negative correlation, R 2 = 0.88, P = 0.01753). 
Deeper analysis among the Proteobacteria revealed that the 
Alphaproteobacteria was the sole class showing significant 
correlations with cadmium measures (R 2 = 0.85; 
P = 0.02718). Moreover, the phylum OP10 exhibited a sig- 
nificant negative correlation with calcium concentrations 
(negative correlation, R 2 = 0.90; P = 0.01395) and the phy- 
lum TM7 was negatively correlated with DOC values (nega- 
tive correlation, R 2 = 0.78; P = 0.04813). 



At the genus level, trends indicated that the intercon- 
nected and disturbed communities of Lar and Bar shared 
two genera with similar co-tolerances to various metalloids: 
Acidisphaera and Acidocella. Furthermore, the community 
of Lar showed additional significant correlations as revealed 
by the genera Acidiphilium, Acidobacteria_GPl Unclassified, 
Kozakia and the genus Rhodovarius. Interestingly, the Lar 
bacterial consortium hosted the genus Tanticharoenia 
which was also found in Tur, the reference community for 
contamination. 

Physicochemical correlation analyses at the structure level 

The results recovered from anova analyses at the structure 
level by assessing the correlations between the relative 
sequences abundance and the physicochemical parameters 
are summarized in Table 7. At the structure level, the cor- 
relations found considering the dominant phyla were com- 
parable to the ones found at the diversity level, with the 
exception that additional correlations were found with 
DOC. Therefore, the relative abundance of sequences 
belonging to the Actinobacteria was negatively influenced 
by the cadmium gradient, while favored by higher dis- 
solved organic carbon content and neutral pH (Cd: nega- 
tive correlation, R 2 = 0.89; P = 0.01605, DOC: R 2 = 0.80; 
P = 0.04009, pH: R 2 = 0.85; P = 0.02522). Within the 
class Actinobacteria, the Actinomycetales still mainly 
explained such a negative correlation with cadmium (nega- 
tive correlation, R 2 = 0.82; P = 0.03407). As found at the 
diversity level, the structure of the phylum Proteobacteria 
significantly correlated with cadmium concentrations, and 
negatively correlated with alkaline pH environments and 
high DOC values (Cd: R 2 = 0.82; P = 0.03325, DOC: neg- 
ative correlation, R 2 = 0.80; P= 0.04116, pH: negative 
correlation, R 2 = 0.86; P = 0.02199). Interestingly, the 
class Alphaproteobacteria did not support the correlation 
found between the structure of the phylum Proteobacteria 
and the cadmium gradient as observed at the diversity 
level. 

At the genus level, Lar and Bar still shared the two same 
genera correlated with multiple metalloids, namely the 
Acidisphaera and the Acidocella. Moreover, Lar community 
still showed significant correlations considering the genera 
Kozakia and Rhodovarius. However, no more correlation 
with cadmium was found for the disturbed communities of 
Lar and Tur. In contrast, the undisturbed community of 
Opa still showed a negative correlation with cadmium con- 
sidering the genus Pelagibacter, a genus also shared with 
the undisturbed community of Das. Moreover, an addi- 
tional negative correlation between cadmium and the rela- 
tive sequences abundance of the genus Beijerinckia was 
found for Opa, and shared with Das and Bar communities. 
The genus Zunongwangia, newly retrieved from the 
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Table 5. Survey of the ten dominant genera found within each community considering the taxonomical structure (based on the relative sequences 
abundance). 
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Lar, Arnoux Lake; Bar, Arnoux Bay; Das, Dasserat Lake; Tur, Turcotte Lake; Opa, Opasatica Lake. 



community of Opa, also exposed such a negative influence 
induced by the cadmium content. 

Interphylum correlation analyses using the relative 
OTUs abundances revealed that the species diversity 
between the phyla Proteobacteria and Bacteroidetes was 
negatively correlated (negative correlation, R 2 = 0.96; 
P = 0.003002), as for the Proteobacteria and the Actino- 
bacteria (negative correlation, R 2 = 0.95; P = 0.00478), 
while a significant positive correlation was detected 
between the OTUs diversity of the Bacteroidetes and the 
Actinobacteria (R 2 = 0.86; P = 0.02321). Considering the 
communities' structure as assessed by the relative 
sequence abundance for each OTU, similar interphylum 
correlations were observed. In this respect, the structure 
of the phyla Proteobacteria and Bacteroidetes were nega- 
tively correlated (negative correlation, R 2 = 0.95; 
P = 0.004256), as for the Proteobacteria and the Actino- 
bacteria (negative correlation, R 2 = 0.96; P = 0.003340), 
while a significant positive correlation was detected 



between the structure of the Bacteroidetes and the Actino- 
bacteria (R 2 = 0.87; P = 0.02152). 

Bacterial network analysis 

Among the 50 dominant species, only 31 exhibited strong 
and significant correlations. Based on these significant cor- 
relations, three abiotic parameter-associated networks and 
six bacteria-exclusive networks were built (see Fig. 7). The 
first abiotic parameter-associated network identifies cad- 
mium, which was found to be positively correlated with the 
OTU Swaminathania, while negatively correlated with the 
abundance of the species Fodinibacter, Thiorhodospira, 
Phenylobacterium and Pelagibacter, Flavobacterium, and 
Arcicella. The second parameter-associated network identi- 
fies the DOC parameter positively correlated with the rela- 
tive abundance of Ferrimicrobium. Ferrimicrobium was 
found to be negatively correlated with the Simiduia abun- 
dance, which in turn was positively correlated with the 
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Table 6. Survey of the significant correlations found considering the OTUs of each community and defining the communities' diversity. 





Al 


Cd 


Cu 


Fe 


Pb 


Zn 


Mn 


Na 


S 


Ca 


Mg 


DOC 


pH 


Temp. 


H' 


NS 


MS 


NS 


NS 


IMS 


IMS 


IMS 


NS 


IMS 


IMS 


IMS 


*(0 90) 


NS 


IMS 


P1 


NS 


t(0 95) 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


*(0.87) 


NS 


P1 .1 


NS 


*(0.88) 


IMS 


NS 


|\|S 


[MS 


IMS 


[MS 


NS 


[MS 


[MS 


IMS 


IMS 


IMS 


P2 


NS 


*(0 90) 


NS 


NS 


IMS 


IMS 


IMS 


IMS 


NS 


IMS 


IMS 


IMS 


*(0.88) 


IMS 


P2.1 


NS 


*(0 85) 


NS 


NS 


IMS 


[MS 


IMS 


NS 


IMS 


IMS 


[MS 


IMS 


NS 


IMS 


P3 


NS 


NS 


NS 


NS 


IMS 


IMS 


IMS 


NS 


IMS 


IMS 


IMS 


*(0 78) 


NS 


IMS 


P4 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


*(0 90) 


NS 


NS 


NS 


NS 


G1 


t(0 93) 


NS 


*(0 91 ) 


+(1 00) 


WO 98) 


*(0 89) 


*(0 81 ) 


[MS 


*(0 81) 


NS 


IMS 


IMS 


IMS 


t(0 93) 


G2 


*(Cl 85) 


NS 


*(0 8?) 


WO 98) 


t(1 00) 


*(0 79) 


IMS 


IMS 


NS 


IMS 


IMS 


IMS 


IMS 


*(0.86) 


G3 


t(0 95) 


NS 


td on) 


*(0 87) 


*(0.80) 


t(0 99) 


WO 93) 


NS 


*(0 91) 


NS 


IMS 


IMS 


*(0 77) 


t(0 98) 


G4 


i(Cl 96) 


NS 


■wo 9?) 


tn oo) 

-t-\ 1 uu i 


WO 96) 


*(0 91 ) 


*(0 85) 


NS 


*(0 87) 


IMS 


IMS 


IMS 


NS 


t(0 94) 


G5 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


*(0 90) 


NS 


NS 


NS 


NS 


G6 


NS 


NS 


NS 


NS 


IMS 


[MS 


NS 


[MS 


NS 


*(0 90) 


[MS 


IMS 


[MS 


IMS 


G7 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


*(0.90) 


NS 


NS 


NS 


NS 


NS 


NS 


G8 


t(0 86) 


[\|S 


*(0 83) 


WO 99) 


+(1 00) 


*(0.80) 


IMS 


NS 


IMS 


IMS 


[MS 


IMS 


NS 


*(0.88) 


G9 


NS 


*(0 82) 


NS 


NS 


IMS 


IMS 


IMS 


NS 


IMS 


IMS 


IMS 


*(0 81 ) 


*(0 90) 


IMS 


G10 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


*(0.89) 


*(0.85) 


NS 


NS 


NS 


G11 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


*(0.89) 


NS 


NS 


NS 


NS 


NS 


NS 


G12 


NS 


*(0.83) 


NS 


NS 


NS 


NS 


NS 


*(0.79) 


NS 


NS 


NS 


NS 


*(0.85) 


NS 


G13 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


*(0.90) 


NS 


NS 


NS 


NS 


G14 


NS 


NS 


* (0.86) 


NS 


NS 


*(0.83) 


NS 


NS 


NS 


NS 


NS 


NS 


*(0.88) 


*(0.85) 


G15 


t(0.95) 


NS 


t(0.94) 


1(0.99) 


t(0.96) 


t(0.92) 


*(0.85) 


NS 


*(0.88) 


NS 


NS 


NS 


NS 


f(0.96) 


G16 


NS 


*(0.82) 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


G17 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


*(0.78) 


NS 


NS 


G18 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


NS 


*(0.78) 


NS 


NS 


NS 


NS 



*P-value between 0.01 and 0.05. 
fP-value between 0.001 and 0.01 . 
JP-value between 0 and 0.001 . 
NS = P-value > 0.05; (R 2 ). 

Data highlighted in grey are negative correlations. The column for potassium was removed since no significant correlations were found. 

Al, Aluminum; Ca, Calcium; Cd, Cadmium; Cu, Copper; Fe, Iron; Mg, Magnesium; Mn, Manganese; Na, Sodium; Pb, Lead; S, Sulfur; Zn, Zinc; DOC, 

Dissolved Organic Carbon; Temp., Temperature; H', Shannon Index. 

Phyla: P1 = Actinobacteria , P1.1 = Actinomycetales, P2 = Proteobacteria , P2.1 = Alphaproteobacteria , P3 = TM7, P4 = OP70; Genera: G1 = Acidi- 
sphaera, G2 = Acidiphilium, G3 = Acidobacteria_GP1 , G4 = Acidocella, G5 = Citreimonas, G6 = Fodinicurvata, G7 = llumatobacter, G8 = Kozakia, 
G9 = Laribacter, G10 = OP10_Undassified, G11 = Orizihumus, G12 = Pelagibacter, G13 = Ponticoccus, G14 = Rhodoferax, G15 = Rhodovarius, 
G16 = Tanticharoenia, G17 = TM7_Undassified, G18 = Wenxinia. 



TM7_unclassified OTU. The third abiotic parameter-associ- 
ated network put in evidence the negative correlation 
between the pH parameter and the abundance of Acidisph- 
aera. Besides these particular abiotic parameter-associated 
networks, networks based exclusively on correlations found 
between bacterial species were also identified. The first 
bacteria-exclusive network regrouped OTUs of the 
Actinobacteria (llumatobacter, Klugiella, Janibacter, and 
Intrasporangium), few members of the Proteobacteria 
(Polaromonas, Oxalicibacterium, and Afifella) and a sole 
representative of the Bacteroidetes (Zunongwagia) . The sec- 
ond bacteria-exclusive network regrouped species from the 
Actinobacteria (Rhodoferax, Pseudorhodoferax, and Hyalan- 
gium) and members of the Bacteroidetes (Fulvivirga and 
Persicobacter) . The third bacteria-exclusive network 
involved only Actinobacteria member interactions, namely 



with OTUs Methylovorus, Beijerinckia, Ideonella, and Der- 
xia. The three last bacteria-exclusive networks consisted in 
single correlations between two species: Seinonella - Roseo- 
monas, OP10_unclassified - Roseomonas, and Parasegetib- 
acter - Catellibacterium. 

Discussion 

Through the identification of taxonomic interactions 
impacted by similar selective pressures, our study aimed to 
shed light on the importance of analyzing communities 
from a phylogenetic perspective. However, we also demon- 
strated that the inclusion of a quantitative approach pro- 
vides even deeper insight into the evolutionary forces that 
shape the dynamics of key taxonomic interaction networks 
in bacterial communities. Our first objective was to assess 
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Table 7. Survey of the significant correlations found considering the sequences of each community and defining the communities' structure. 
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whether natural selection drove parallel adaptation at the 
community taxonomical diversity and structure levels on 
two geographically isolated bacterial communities under 
similar selective pressure. The second objective was to test 
whether a polymetallic gradient exposure in three intercon- 
nected lacustrine bacterial communities triggered signifi- 
cant taxonomical membership and abundance shifts. 
Finally, the third objective was to identify the related envi- 
ronmental factor driving community changes in the stud- 
ied system. 



Despite recent progress in characterizing the diversity of 
freshwater microorganisms, drivers of changes in commu- 
nity composition have not yet been well characterized. 
Anthropogenic environmental disturbance may have a 
great impact on the functional diversity of bacterial pro- 
cesses. Given such disturbances are increasingly common, 
there is a need for accurate predictive models regarding 
their impact on bacterial community ecology and taxo- 
monic diversity (Konopka 2009). More specifically, as taxo- 
nomic diversity is suspected to influence the metabolic 
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Figure 7 Bacterial networks of the 50 dominant species and the physicochemical parameters of interest (Cd = cadmium; DOC = Dissolved organic 
carbon; pH = pH measure). The nodes represents the OTUs with 97% identity at a P-value of <0.01 and each node color indicates to which phylum 
the OTU belongs (yellow = Actinobacteria, blue = Bacteroidetes, purple = Proteobacteria, green = Firmicutes, orange = OP10, and grey = 7M7). 
Each strongly significant correlation is an edge represented by a grey line. Edges indicative of a negative significant correlation are represented by a 
red line. 



pathway interacting networks, a first step is to identify key 
taxonomic interaction networks that are lost in the face of 
strong environmental perturbation. Using high-throughput 
sequencing of 16S ribosomal subunit gene libraries, this 
study successfully identified environmental parameters 
with the power to drive key taxonomic interaction net- 
works in bacterial communities, both at the community 
membership (diversity) and composition (structure) levels. 
In doing so, we characterized the impact of a polymetallic 
contamination gradient over a system of both independent 
and interconnected lacustrine bacterial communities. 

Among the collected data, geochemical analysis of lake 
water samples indicated that the heavy metal content in 
Opa was similar to the one found in Das, while the heavy 
metal concentrations in Tur were similar to those measured 
in Bar. Consistent with the geochemical data, PCoA, Venn 
and phylogenetic analyses at the microbial diversity level 
supported the clustering of undisturbed communities on 
one side and disturbed communities on the other. In this 
respect, the clustering of the geographically isolated com- 
munity of Tur with the interconnected communities of Bar 
and Lar revealed that parallel adaptation can occur at the 
taxonomic level between independent communities if 
impacted by a similar selective pressure. Such a parallel 
change of interacting taxonomical networks suggests that 



community adaptation is the net effect of individual mem- 
bers' successful and unsuccessful acclimation/adaptation to 
environmental perturbation, as proposed by DeAngelis 
et al. (2010). Therefore, parallel changes of interacting 
taxonomical networks would suggest in turn that HGT of 
loci corresponding to heavy metal resistance were, at the 
most, limited to closely related species, as argued by An- 
dam and Gogarten (2011). Thus, our results suggest that 
two independent bacterial communities living in similar 
environmental conditions exhibit the same overall taxo- 
nomical membership (i.e. members in both communities 
belonging to the same high taxonomic ranks). 

At the structure level (i.e. OTU relative abundance), 
PCoA, Venn and phylogenetic analyses supported the clus- 
tering of the communities into three groups according to 
the cadmium gradient, thus exerting the fact that cadmium 
was the most lethal heavy metal measured in this study. As 
such, in relation to relative taxon abundance, weighted 
structural analyses grouped communities Tur and Bar into 
an intermediate metallic contaminated communities clus- 
ter and further demonstrated that higher cadmium concen- 
trations thoroughly disturbed Lar community structure 
(Fig. 4B). This result was striking because Lar and Bar lakes 
are connected. Overall, these results first bring to light evi- 
dence that community differentiation at the taxonomical 
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diversity level reflects the presence/absence of a perturba- 
tion, while the community structure level accurately reflects 
the magnitude of the perturbation. 

Our study has therefore enlightened how parallel com- 
munity adaptation resulting from a similar environmental 
perturbation may have occurred within the original core 
microbiome due to selection of a heavy metal resistant core 
OTUs combined to the loss of most sensitive bacterial spe- 
cies. Consequently, all the OTUs shared across the five lakes 
may be considered as part of the original core microbiome. 
Then, the OTUs shared between Opa and Das may be iden- 
tified as the healthy core microbiome mirroring the absence 
of contamination while the persistent species shared 
between the polluted sites make up the heavy metal resistant 
core microbiome selected from the original core microbiome. 
A comparable bacterial community adaptation pattern was 
observed for copper contamination, where membership 
shifts translated by an increase of Proteobacteria (Turpeinen 
et al. 2004). In this study, a similar trend was observed: 
Proteobacteria were widely prevalent among the disturbed 
communities, accounting for more than 84% of the OTUs 
and their relative abundance reached more than 92% of the 
sequences, while their proportion dropped to less than 
60% of the OTUs and their abundance was less than 71% 
of the sequences in the undisturbed communities. As 
revealed by the correlation analyses, the diversity among 
the Proteobacteria was mainly favored by higher cadmium 
level and acidic pH values. The Proteobacteria dominance 
was also influenced by these two previous factors, and by 
the dissolved organic content. Such finding could be 
explained by either the presence of heavy metals resistant 
species among the Proteobacteria and/or the fact that this 
particular phylum encompasses more generalist species, 
thus conferring competitive traits to the whole phylum in 
case of a perturbation (Barberan et al. 2012). Furthermore, 
the Alphaproteobacteria OTUs diversity was found to 
underlie such a correlation between Proteobacteria and cad- 
mium. Interestingly, among the Alphaproteobacteria, sev- 
eral genera showed positive correlations with Fe, the 
highest metalloid measured in the disturbed sites, but also 
the most biologically important metal cation owing to its 
importance in the electron transport chain (Nies 1999). 
Furthermore, numerous Alphaproteobacteria genera exhib- 
ited positive correlations with Al, Cu, Pb, Zn, and Mn. This 
observation illustrates Alphaproteobacterial resiliance to 
polymetallic contamination and may account for the high 
Chao richness and the number of OTUs observed (S Q b s ) in 
Lar. Alphaproteobacterial genera encompassed Acidiphili- 
um, Acidimicrobium, Acidisphaera, Acidocella, Kozakia, and 
Rhodovarius, which were found to be favored at both the 
diversity and structure levels in Lar. Interestingly, the Al- 
phaproteobacteria are known to be highly successful aquatic 
lineages. Competitive at low nutrient levels and resistant to 
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grazing by organisms of higher trophic levels (Newton 
et al. 2011), they also bear genomes enriched with genes 
belonging to xenobiotic degradation pathways (Debroas 
et al. 2009). Moreover, heterotrophic acidophiles belonging 
to the Alphaproteobacteria class such as the genera Acidiphi- 
lium, Acidocella, and Acidomicrobium appeared to be widely 
distributed in metal-rich acidic environments as they cata- 
lyze processes implied in the iron cycling (Harrison 1984; 
Johnson and McGinness 1991; Kiisel et al. 1999, 2002; 
Hallberg and Johnson 2001, 2003). Our results therefore 
confirm the overall opportunistic ability of the Proteobacte- 
ria, particularly the class of Alphaproteobacteria, commonly 
identified in lake epilimnion. 

The phylum Actinobacteria is also highly abundant in 
freshwater environments (Newton et al. 2007, 2011). In 
this study, the phylum Actinobacteria was severely impaired 
both at the diversity and structure levels, by exposure to 
cadmium, acidic pH, and low DOC content, especially con- 
sidering the Actinomycetales order. Such results were 
expected as cadmium is known to be among the most 
harmful heavy metal for microorganisms (Nies 1999). 
Decrease in pH is known to induce a decline in the bacte- 
rial structural and functional diversity (Anderson et al. 
2009). In a detailed analysis of tribes within the phylum Ac- 
tinobacteria, Newton et al. (2007) recorded distribution 
patterns along pH gradients revealing that lake pH is a 
strong predictor of the community composition in this 
particular phylum. Interestingly, DOC is a parameter 
directly related to the organic matter decomposition activ- 
ity of the bacterial community. In this respect, Actinobacte- 
ria are believed to play a major role in organic matter 
degradation (Steger et al. 2007). Taken together, it suggests 
that the significant correlation between the low values of 
DOC and the drastic decrease of actinobacterial strains 
may be attributed to the loss of related metabolic pathways 
induced by heavy metal contamination. 

Interphylum correlation analyses revealed that the Ac- 
tinobacteria were shared an inverse relationship with the 
Proteobacteria. An inverse correlation was also detected 
between the Bacteroidetes and the Proteobacteria, despite 
the fact that the Bacteroidetes showed only marginally sig- 
nificant correlations with both cadmium and pH when 
considering the structure and diversity levels (data not 
shown). We can therefore conclude that there were strong 
patterns of co-variation between phyla, suggesting that in- 
terphylum interactions are affected in parallel across com- 
munities when impacted by a similar perturbation. This 
result suggests in turn that taxonomic interphylum interac- 
tions are likely conserved. 

The significant correlations found with the anova analy- 
ses between environmental parameters and the phylogenic 
survey strongly suggest that both cadmium and pH exert 
the widest influence over the diversity and structure of 
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bacterial communities. Furthermore, network analyses 
highlight that the pH parameter was negatively correlated 
with the Acidisphaera, meaning that acidic waters favor 
this particular genus, consistent with its extensive isola- 
tion from acidic environments (Hiraishi et al. 2000; John- 
son et al. 2001). However, as revealed by the network 
analyses, cadmium is by far the parameter that exerted 
the widest negative impact over the abundance of the 
dominant species. Effectively, multiple negative correla- 
tions were found with cadmium, negatively affecting the 
abundance of 12% of the dominant species. Such findings 
strongly suggest that cadmium disturbs a significant pro- 
portion of the microbial biosphere as hosted species inter- 
act through taxonomic networks. Because dominant 
species have a large impact on the cycling of important 
nutrients (Malmstrom et al. 2004; Campbell et al. 2011), 
our results strongly suggest that cadmium would severely 
impede a variety of ecosystem functions. As a conse- 
quence, highly impacted sites have lower DOC content 
due to the reduced bacterial decomposition activity (Do- 
elman and Haanstra 1979; Kelly and Tate 1998; Kelly 
et al. 1999). Such observations may explain that the genus 
Ferrimocrobium, although commonly found in acid mine 
drainages due to its capacity for ferrous iron oxidation, 
showed a lower abundance in the highly impacted sites 
where the lowest DOC values would result from weaken 
bacterial activities. This particular genus may therefore 
regroup sensitive bacterial species with great potential for 
assessing the intensity of a disturbance over the bacterial 
activities sustaining ecosystem services such as the organic 
matter decomposition. 

Overall, our study brought evidence confirming the 
hypothesis that following a strong perturbation, sensitive 
bacterial lineages are eliminated, thus allowing opportunis- 
tic species to dominate the emptied niches. The disappear- 
ance of sensitive species gradually drives the ecosystem to 
an overall decrease in biodiversity, which in turn favors 
biotic homogenization (McKinney and Lockwood 1999). 
To this respect, detecting dynamic changes of both con- 
served taxonomic interphylum interactions and species 
level interaction networks provides new insights to under- 
stand further the consequences of a homogenization phe- 
nomenon that have been widely mentioned in the literature 
(Kandeler et al. 1996; Pennanen et al. 1996; Aoyama and 
Nagumo 1997; Ciller et al. 1998; Kelly et al. 1999). 

In conclusion, our results demonstrate that heavy metal 
contamination caused significant shifts within lacustrine 
bacterial communities both at the membership and struc- 
ture levels, and drove parallel adaptation of two indepen- 
dent bacterial communities impacted by a similar selective 
pressure. Determining the extent to which ecologically rel- 
evant traits are phylogenetically conserved in such a paral- 
lel way is critical to understand the parameters that 



control the community phylogenetic structure, and their 
evolutionary consequences. In this respect, the parallel 
changes of interacting taxonomical networks would sug- 
gest that HGT of heavy metal resistance genes were limited 
to closely related species, as argued by Andam and Gogar- 
ten (2011). Otherwise, different patterns of interacting 
taxa would have been observed among independent con- 
taminated lake communities. Therefore, investigating the 
functional repertories that were favored (gained) or disfa- 
vored (loss) in contaminated lakes will allow assessing fur- 
ther whether the taxonomic diversity of bacterioplankton 
communities primarily influenced the pathways of meta- 
bolic response to different environmental factors, as sug- 
gested by Comte and del Giorgio (2011). This study is 
underway. 

Finally, cadmium and acidic pH were identified as the 
physicochemical parameters that exerted the strongest 
influence on the whole community taxonomic network. 
These two physicochemical parameters favored the Alpha- 
proteobacteria, which were observed to have a functional 
repertory enriched in detoxifying genes (Debroas et al. 
2009). Consequently, with respect to pollution control, this 
phylum may regroup species with valuable properties with 
regard to (i) hastening the mobilization of metals (with a 
possible application in biomining) (Vails and de Lorenzo 
2002), (ii) designing metal- tolerant strains that are better 
adapted to perform biodegradation of organic pollutants, 
and (iii) metal bioremediation or mitigation through the 
breeding of natural or engineered strains. 

Such insights therefore clearly demonstrate the need to 
study bacterial community diversity and structural changes 
triggered by abiotic contaminant exposure from an ecologi- 
cal and evolutionary point of view. Heavy metals affect the 
bacterial communities that occupy different ecological 
niches and participate in a variety of biogeochemical pro- 
cesses and ecosystem functions. As the world's freshwater 
resources are under severe stress, next-generation sequenc- 
ing and bioinformatics are promising avenues to explore 
the functional interacting network adaptability of bacterial 
communities, and in turn, predict their capacity to main- 
tain ecosystem homeostasis when facing rapid environmen- 
tal changes. 
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